addpath(genpath(pwd))
clear
%% Model specific part
unitFree = 1;
[f,x,xp,y,yp,symparams,Phi,h1,nomSDF_cup,heteroShock,r_cu] = DSGEmodel(unitFree);

% The number of endogenous states (mx) and number of lagged controls (myx)
mx  = 0;
myx = 1;
ne  = 5;
%% The analytical derivatives of the model
% The dimension of the model
ny  = length(y);
nx  = length(x);

% The parameters of the model
modelParams = cell(1,length(symparams));
for i=1:length(symparams)
    modelParams{i} = char(symparams(i));
end

% Remove exogenous shocks
fall         = f;    
f            = f(1:ny+mx,1);
%% Generating 'getModelDeriv.m'
% Folder were the file is stored
positionFiles     = [pwd,'\ModelSpecification'];

% Function for computing the residuals for projection 
strucOfArrays = struct('f',f);
nx1           = mx + myx;
nameOfFunction = 'ResidualEquationProj_raw.m';
generateMfuncResidualEquation_file(modelParams,x,y,xp,yp,nx1,ne,h1,heteroShock,strucOfArrays,positionFiles,nameOfFunction)
DispSymMatrixMatlab_file_vector(r_cu,'r_cu','r_cu_display.txt')
DispSymMatrixMatlab_file_vector(nomSDF_cup,'nomSDF_cup','nomSDF_cup_display.txt')
% then manually modify ResidualEquationProj_raw.m to get ResidualEquationProj.m

% For computing bond prices
generateMfuncSaveStatesControls_file(modelParams,x,y,xp,yp,nx1,ne,h1,heteroShock,positionFiles)
nameOfFunction = 'ResidualEquationProjBondPrices.m';
BondEq = [char(nomSDF_cup),'*exp(p1_cup)'];
generateMfuncResidualEquationBondPrices_file(modelParams,x,y,xp,yp,nx1,ne,h1,heteroShock,BondEq,positionFiles,nameOfFunction)

% For computing real bond prices
nameOfFunction = 'ResidualEquationProjBondPricesReal.m';
syms pai_cup;
SDF_cup = simplify(nomSDF_cup*exp(pai_cup));
BondEq = [char(SDF_cup),'*exp(p1_cup)'];
generateMfuncResidualEquationBondPrices_file(modelParams,x,y,xp,yp,nx1,ne,h1,heteroShock,BondEq,positionFiles,nameOfFunction)

% For computing expected values of controls
nameOfFunction = 'ResidualEquationProjExpectedControl.m';
generateMfuncResidualExpectedControl_file(modelParams,x,y,xp,yp,nx1,ne,h1,heteroShock,positionFiles,nameOfFunction)

% For computing equity prices
%nameOfFunction = 'ResidualEquationProjEquityPrice.m';
%syms muz_cup
%discountFunc   = simplify(nomSDF_cup*exp(pai_cup)*exp(muz_cup));    %adding exp(muz_cup) to account for growth
%generateMfuncResidualEquityPrice_file(modelParams,x,y,xp,yp,nx1,ne,h1,heteroShock,discountFunc,positionFiles,nameOfFunction)
